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[57] ABSTRACT 

An electrophysiological activity estimation method is used 
to estimate active areas of a selected part of a living body 
(e.g., human brain). Based on electromagnetic field distri- 
bution measured on a surface of the selected part of the 
living body as well as its shape, the method assumes each 
equivalent current dipole as a source of electromagnetic field 
distribution. At first, grid points are created to cover the 
selected part of the living body (e.g., head). Using the 
pseudoinverse, the method performs estimation of dipole 
distribution for equivalent current dipoles. Then, the equiva- 
lent current dipoles are subjected to sorting based on a 
priority order. An evaluation function (e.g., structural risk) is 
calculated in accordance with a square error between the 
measured electromagnetic field distribution and a theoretical 
value of electromagnetic field distribution caused by the 
equivalent current dipole, as well as a number of measuring 
points and a number of equivalent current dipoles. Estima- 
tion and calculations are repeated while changing a number 
of equivalent current dipoles in accordance with the priority 
order. Then, noise distribution is estimated in response to a 
minimum value of the evaluation function. So, the dipole 
distribution is estimated using the noise distribution. Thus, 
it is possible to obtain estimation results which are repre- 
sented by a direction and magnitude of the equivalent 
current dipole as well as an appropriate number of equiva- 
lent current dipoles. 

19 Claims, 5 Drawing Sheets 
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ELECTROPHYSIOLOGICAL ACTIVITY 
ESTIMATION METHOD 

BACKGROUND OF THE INVENTION 

1. Field of the Invention 5 
This invention relates to electrophysiological activity 

estimation methods which estimate equivalent current 
dipoles (i.e., ECD) based on the electromagnetic field dis- 
tribution which is measured on a surface of a living body 10 
such as a scalp of a human. This invention is based on patent 
application No. Hei 9-19619 filed in Japan, the content of 
which is incorporated herein by reference. 

2. Prior Art 

Conventionally, the equivalent current dipole estimation 15 
method is known as the method to estimate active areas of 
a human brain on the basis of electromagnetic field distri- 
bution measured on a scalp of a human. By observing active 
areas of the human brain, it is possible to obtain information 
with regard to high-order functions of the brain as well as 20 
diseased parts of the brain. 

Now, a general method to calculate an equivalent current 
dipole will be described with respect to estimation of active 
areas of the brain. When a nerve cell (or neuron) of the brain 
is excited by an impulse (or stimulus) given thereto from an 25 
external device, pulse-like current flows across a connection 
(or axon) connecting neurons together. Due to occurrence of 
the current, weak electromagnetic field is caused to occur 
around a scalp of the human. A source of the electromagnetic 
field is subjected to modeling using current elements called 30 
equivalent current dipoles (simply called "dipoles"). By 
estimating six parameters of the dipole, it is possible to 
estimate an active area of the human brain. Herein, the six 
parameters contain three parameters for designation of an 
area, two parameters for designation of a direction and one 35 
parameter for representation of magnitude (or intensity) with 
respect to the dipole. 

Two methods are known as the method to estimate the 
parameters of the dipole. Herein, the six parameters of the ^ 
dipole are unknown while a measured value of the electro- 
magnetic field is represented by an equation as follows: 



y-tyxv-b > > yn) T 



[Equation 1] 



A theoretical value is calculated using the dipole model 45 
represented by an equation as follows: 

/(6/Mf i(»j>. Ufy, • • • , UQjV. [Equation 2] 

So, a square residual is calculated as follows: 50 



measured value and theoretical value of the electromagnetic 
field is represented as follows: 

y,f 

Incidentally, an expression of (. . .) T represents the transpo- 
sition. Further, a coordinate of the measuring point i is 
represented by 



while a dipole parameter at the measuring point j is repre- 
sented by an equation as follows: 



[Equation 5] 



In addition, the position, direction and magnitude of the 
dipole at the measuring point j are respectively represented 
by symbols as follows: 

As the method to estimate the foregoing dipole 
parameters, it is possible to use the nonlinear optimization 
algorithms such as the Levenberg-Marquardt method and 
simplex method. Details of the nonlinear optimization algo- 
rithms are shown by "paper 1" which is a book entitled 
"Nonlinear Programming" which is written by Mr. Hiroshi 
Yamashita and Mr. Hiroshi Konno and is provided by 
Science Technology Association of Japan on 1987, for 
example. Further, the method using the pseudoinverse of the 
matrix (or Moore-Penrose inverse) is known as an example 
of the other methods for the current dipole estimation, 
especially as an example of the estimation method for 
"distributed" active areas of the living body. According to 
this method, the theoretical equation of the electromagnetic 
field is represented as follows: 



[Equation 6] 



So, it is possible to use the characteristic of the above 
equation that elements regarding the magnitude of the dipole 
are coordinated by the linear function. Concretely speaking, 
the shape of the brain is approximated by the set of the 
polyhedrons while the position and direction of the dipole 
are fixedly located on the polyhedron, so that the magnitude 
q, of the dipole is calculated as an unknown parameter. In 
this case, the function of "F* is represented by the matrix, so 
the aforementioned equation 6 is rewritten in the matrix 
form as follows: 



n 



[Equation 3] 



The aforementioned method estimates the parameter mini- 
mizing the square residual, which is represented by an 
equation as follows: 



[Equation 4] 



In the aforementioned equations, "ri" denotes a number of 
measuring points; "m" denotes a number of dipoles; "y/* 
denotes a measured value of the electromagnetic field mea- 
sured at a measuring point i; u f" denotes a theoretical value 
of the electromagnetic field at the measuring point i. In 
addition, a vector corresponding to combination of the 



/=Fq 



[Equation 7] 



For example, the electroencephalograph (i.e., EEG which 
55 represents a brain wave measuring instrument) is used to 
measure the electric field distribution (or electric potential 
distribution), based on which the dipole is estimated. In this 
case, calculation for F is made with regard to elements i, j 
in accordance with an equation as follows: 

60 



y 4,rZj|i 



[Equation 8] 



65 In the above equation, "Y" denotes the spherical harmonics; 
symbols r,0,<J> show the position of the measuring point in 
polar coordinates while symbols r'^',^ 1 show the position of 
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the dipole in polar coordinates. A symbol V denotes a 
differential operator while a mathematical expression of 
R(r,r') corresponds to the function regarding r and r', which 
is determined by the boundary condition. Using a SQUID 
(an abbreviation for "Superconducting Quantum Interfer- 
ence Device"), it is possible to measure the magnetic field 
distribution, based on which the dipole is estimated in 
accordance with an equation as follows: 

f - Po n; ' e > X fo {Equation 9] 

where ^ is a magnetic permeability. In addition, n ( . denotes 
a normal vector with respect to a surface of a head model at 
a measuring point i. A square error is represented by an 
equation as follows: 

E=\\y-F4 2 [Equation 10] 

Herein, a vector q which minimizes the above square error 
is represented, using a pseudoinverse F* of F, by an equation 
as follows: 

q=F*y [Equation 11] 

The dipole estimation method using the pseudoinverse is 
disclosed by "paper 2" corresponding to Japanese Patent 
Laid-Open Publication No. 6-343613 as well as "paper 3" 
corresponding to Japanese Patent Laid-Open Publication 
No. 6-343614, for example. 

Now, FIG. 6 is a flowchart showing an example of 
processing which the conventional system executes in accor- 
dance with the conventional method using the pseudoin- 
verse. First, a group of grid points are arranged uniformly 
with respect to a diagnosing area, i.e., an area on which 
active areas are estimated in step 51. In step 52, a current 
source is calculated with regard to each of the grid points in 
accordance with the minimum-norm least-squares method. 
In step 53, a group of grid points are moved in proximity to 
a grid point on which a current source having a large value 
exists. Then, a minimum space between grid points is 
selected from among spaces between the grid points. In step 
54, a decision is made as to whether the minimum space is 
less than a prescribed value or not. If the minimum space is 
less than the prescribed value, the system completes the 
proceeding of FIG. 6 because the convergence is obtained. 
If not, the system proceeds back to step 51. 

According to the aforementioned method, the system 
performs estimation while changing the distance between 
the grid points in order to improve a precision of estimation. 
However, the conventional system does not perform esti- 
mation with respect to a number of dipoles. In addition, the 
conventional system is not designed to cope with problems 
due to the noise. 

Next, a description will be given with respect to problems 
which the conventional dipole estimation method suffers 
from. 

Normally, the method using the nonlinear optimization 
algorithm uses one to three dipoles to perform estimation of 
an active area. However, it is impossible to know a number 
of the active areas in the brain in advance. For this reason, 
there are provided dipole models whose number of dipole(s) 
is one, two and three respectively, for example. Using the 
dipole models, the system performs the estimation. So, the 
system employs the best result of the estimation. In this case, 
a square error E can be used as an element of evaluation 
criterion for the quality of estimation. However, there is a 
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tendency that the square error becomes smaller as a number 
of dipoles is increased more. For this reason, there is a 
problem that a number of dipoles cannot be detected with 
accuracy. FIG. 7 shows an example of variations of the 

5 square error of the estimation which is made with respect to 
the magnetic field distribution produced by two dipoles 
while changing a number of model dipole(s) from one to 
eleven. In FIG. 7, the vertical axis corresponds to the 
logarithmic notation of the mean square error (or mean 

10 square residual, i.e., MSR) while the horizontal axis repre- 
sents a number of model dipoles. According to the graph of 
FIG. 7, it is observed that the square error is decreased as a 
number of model dipoles is increased. So, if the square error 
is used as the evaluation criterion, although a true number of 

IS dipoles is two, the best result of the estimation is obtained 
with respect to the case of eleven model dipoles. This is a 
wrong result. In addition, the method using the nonlinear 
optimization algorithm requires a long time to perform the 
estimation one time. Further, it is necessary to perform the 

20 estimation several times while changing the number of the 
dipoles. Namely, the conventional method requires a large 
load to calculations as well as a long time for the estimation. 

Another conventional method using the pseudoinverse 
arranges several hundreds to 10,000 dipoles with regard to 

25 each element of the polyhedrons approximating the shape of 
the brain, wherein magnitude of each dipole is calculated. 
According to this method, the estimation can be made using 
multiplication of matrix only. So, this method requires a 
short time for the estimation. Even if a neural activity lies in 

30 an extended area, the method is capable of providing infor- 
mation about the size of the area. The method using the 
pseudoinverse is designed to use the linear model. For this 
reason, this method has a drawback corresponding to weak- 
ness against the noise. 

35 Now, a description will be given with respect to reasons 
why the method using the pseudoinverse has such a draw- 
back. 

Using the singular value decomposition (i.e., SVD), the 
pseudoinverse F^ is written by an equation as follows: 

40 

r = V /\U T =V -v,«f [Equation 12J 

kix> Xi 

45 where u^v,- represent column vectors i for matrices U, V; hi 
represents a singular value of the matrix F, so sum is 
performed with respect to singular values which are not 
zero. If measurement data do not contain noise, correct 
dipole distribution can be obtained by multiplication of the 

50 pseudoinverse and measurement data in accordance with an 
equation as follows: 

q = F+y = YtTi V ^-y) [Equation 13] 



where q corresponds to a solution vector having a minimum 
norm which minimizes the square error. However, in the 
case where the measurement data contain noise z, the 
60 equation 13 should be rewritten by using an equation as 
follows: 

y°y^ z [Equation 14] 

65 where y 0 corresponds to measurement data containing no 
noise. So, using the equation 14, the equation 13 is rewritten 
as follows: 



04/22/2004, EAST Version: 1.4.1 



6,073,040 

6 

Iq addition, by using the noise distribution for the esti- 

i - V i- 7 [Equation 15] mation of the dipole distribution, this invention is capable of 

9- y- tyo z - 2-j Xi v * u ; w constructing an architecture of processing for the dipole 

estimation which is robust against the noise. 

5 BRIEF DESCRIPTION OF THE DRAWINGS 

These and other objects of the subject invention will 

In the above equation, a second term of the right side shows become more fully apparent as the following description is 

an extremely large value with respect to a small singular read m ^g 01 of attached drawings wherein: 

value >u. So, the divergence occurs on an estimated value of FIG. 1 is a flowchart showing procedures of an electro- 

the magnitude of the dipole. physiological activity estimation method in accordance with 

So, the problems of the conventional methods can be an embodiment of the invention; 

summarized as follows: FIG. 2 shows an example of a shape of a human head 

In the case of the method using the nonlinear optimization i5 which is approximated using polyhedrons; 

algorithm where the square error is used as the evaluation FIG 3 is a graph snowing a relationship between a value 

criterion for the estimation, it is impossible to estimate a of m evalualion functioo m6 a number of model 

number of dipoles. In the case of the method using the _ _ . , , , r . . . 

pseudoinven* for estimation of the distributed active area, 4 ». a g^ph showing an example of simulation data 

f. r 4U r*u corresponding to a distributed dipole model; 

the divergence occurs on estimation of the magnitude 01 the i- -& v » 

dipole due to the effect of the noise. 20 mG - 5A ^ a S ra P h showing results of estimation which is 

performed based on this invention with respect to the 

SUMMARY OF THE INVENTION simulation data of FIG. 4; 

It is an object of the invention to provide an electrophysi- FIG. 5B is a graph showing results of estimation which is 

ological activity estimation method for estimation of active performed based on the conventional method using the 

areas of a living body (or organism) which is capable of 25 pseudoinverse with respect to the simulation data of FIG. 4; 

estimating a number of dipoles with accuracy and which is FIG. 6 is a flowchart showing procedures of the conven- 

capable of performing noise-robust estimation of the tional dipole estimation method using the pseudoinverse; 

dipoles. and 

According to an aspect of this invention, an electrophysi- 30 FIG. 7 is a graph showing a relationship between a 
ological activity estimation method is applied to estimation number of model dipoles and a square error with respect to 
of active areas of a selected part of a living body (e.g., the conventional dipole estimation method using the non- 
human brain). Herein, each current dipole is assumed as a linear optimization algorithm, 
source of electromagnetic field distribution on the basis of 

measured electromagnetic field distribution measured on the 35 DESCRIPTION OF THE PREFERRED 

human scalp as well as its shape. EMBODIMENT 

At first, grid points are created to cover the selected part Before describing the preferred embodiment of this 

of the living body (e.g., head). Using the pseudoinverse, the invention, a description will be given with respect to a 

method performs estimation of dipole distribution for principle of the invention. 

equivalent current dipoles. Then, the equivalent current 40 This invention is characterized by introducing a new 

dipoles are subjected to sorting based on a priority order. An evaluation function "L" which is used instead of the square 

evaluation function (e.g., structural risk) is calculated in error used by the conventional methods. Herein, the evalu- 

accordance with a square error between the measured elec- a tion function L is established in consideration of the 

tromagnetic field distribution and a theoretical value of complexity of the models. According to an example of the 

electromagnetic field distribution caused by the equivalent 45 evaluation function, for example, the complexity of the 

current dipole, as well as a number of measuring points and model is evaluated using a number of parameters (i.e., a 

a number of equivalent current dipoles. number of dipoles), wherein an excessive use of parameters 

The above estimation and calculations are repeated while is prohibited. Concretely speaking, a structural risk is used 

changing a number of equivalent current dipoles in accor- as an evaluation function L SRM . The structural risk is defined 

dance with the priority order. Then, noise distribution is 50 by equations as follows: 
estimated in response to a minimum value of the evaluation 

function. So, the dipole distribution is estimated using the £ [Equation 16] 

noise distribution. Thus, it is possible to obtain estimation LsRM ~ c sm 
results which are represented by the direction and magnitude 

of the equivalent current dipole as well as an appropriate 55 

number of equivalent current dipoles. c SRM = 1 -c 

By changing the number of equivalent current dipoles as 
a repetition of the estimation and calculations, it is possible 

to change a group of equivalent current dipoles used for the where E denotes a square error between a measured value 

estimation. In other words, after performing the estimation 60 and a theoretical value of electromagnetic field distribution; 

using the pseudoinverse, this invention excludes unneces- n denotes a number of measuring points while m denotes a 

sary dipoles which are regarded as "fake" neural activities. number of dipole parameters (i.e., a number of dipoles). In 

Therefore, as compared with the conventional methods, the addition, t] is a value regarding a probability, so the above 

electrophysiological activity estimation method of this equation is established with a probability of I-tj. Further, 

invention is advantageous that equivalent current dipole 65 d^m) is called a Vapnik-Cbervonenkis (VC) dimension 

estimation can be performed with a small amount of calcu- which corresponds to the quantity used as an index of the 

lations and at a high speed. diversity or expressiveness of the model. Therefore, the VC 
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dimension becomes large together with a number of param- cretely speaking, the sum of the aforementioned equation 13 

eters of the model. In addition, 0 SRM which is used as a is subjected to limitation that the square error is equivalent 

denominator of the above equation becomes small, so a to the noise distribution. In other words, the equation 13 is 

value of L SJiM as a whole becomes large. rewritten as follows: 

This invention is designed to determine that an optimum 5 

model corresponds to a dipole model which minimizes the _^ _y 1 r [Equation 17] 

new evaluation criterion of L^^. The conventional method q ~ y ~ \i V ' y 
suffers from a problem that a true number of dipoles cannot 
be estimated because a square error becomes small as a 

number of parameters of the model becomes large. In 10 where ||z|| is estimated magnitude of the noise distribution, 

contrast, this invention uses the evaluation function L^, By changing the content of the sum as described above, this 

according to which becomes small together with a invention is capable of avoiding over-fitting to the noisy 

number of dipole parameters. So, an excessive use of data > so it is possible to realize the noise-robust estimation, 

parameters brings a disadvantage in an aspect of minimiza- m this case, singular value decomposition is performed on 

tion of ^srm- According to the mechanism of the invention 15 matrices which are determined by position information of 

described above, it is possible to estimate a number of equivalent current dipoles as well as arrangement of mea- 

dipotes. In the new evaluation function h SRM can be suring points of electromagnetic field distribution. Thus, 

assumed as a penalty term which acts to prohibit the using me noise distribution as well as error Mormation, it is 

excessive use of parameters. Incidentally, details of the possible to perform estimation with respect to directions and 

theoretical basis that the evaluation function L^r^ is effec- 20 magnitude of the equivalent current dipoles, for example, 

live for estimation of a number of dipoles is provided by Incidentally, there are two methods for determination of 

"paper 4" which corresponds to V. N. Vapnik, "The Nature foe aforementioned priority order, as follows: 

of Statistical Learning Theory", Springer, N.Y., 1995. i) The priority order is determined to coincide with an 

Other than the aforementioned evaluation function repre- order of absolute values of the equivalent current 

sented by the equation 16, this invention is capable of 25 dipoles. 

employing an evaluation function based on a description ii) An error is detected between measured electromagnetic 

length L MZ>L , for example. In addition, it is possible to use field distribution and electromagnetic field distribution 

the Akaike information criterion as the evaluation which is caused by each equivalent current dipole. So, 

function. Incidentally, the evaluation function based on the priority order is determined to coincide with an 

L MnL can be assumed as a simplified form of the aforemen- 30 order of errors. 

tioned evaluation function of the equation 16. Next, a description will be given with respect to the 

If the conventional methods, especially the method using preferred embodiment of the present invention. FIG. 1 is a 

the nonlinear optimization algorithm, use the equation 16 flowchart showing concrete procedures to execute an elec- 

which allows estimation of a number of dipoles as the trophysiological activity estimation method in accordance 

evaluation function, it is necessary to perform estimation 35 with one embodiment of the invention. Incidentally, the 

many times while changing the number of model dipoles, embodiment refers to the method to estimate active areas of 

and it is necessary to perform comparison of an evaluated a brain on the basis of electromagnetic field distribution 

value of Ljra/ with respect to each estimation. For this measured on a scalp of a living body called "subject", 

reason, there is a problem that much time is required for the However, this invention is not limited to such estimation of 

estimation. In contrast to the conventional methods, this 40 the active areas of the brain. For example, the method 

invention is designed to perform estimation procedures as similar to the method of this invention can be employed for 

follows: e^tiinalioa^^aiv^areas^^^^ 

At first, estimation is performed using the pseudoinverse ejictiio^ 

with respect to dipole distribution. Then, estimation is breast -of ; the r subjects In addition, this invention can be 

performed with respect to a number of dipoles while delet- 45 applied to estimation of active areas of the digestive organ 

ing unnecessary dipoles from the dipole distribution or and muscle as well. 

adding dipoles having a high priority order to the dipole In step 1 of the flowchart of FIG. 1, a loop counter 

distribution. The method using the pseudoinverse is capable variable k is subjected to initial setting at '0*. In addition, the 

of performing dipole estimation by using multiplication of system of this invention loads data regarding electromag- 

matrices only. So, as compared with the method using the 50 netic field distribution measured at measuring points on 

nonlinear optimization algorithm which requires iterative a head of a subject Further, the system loads data regarding 

calculations, the method using the pseudoinverse is capable a shape of the head of the subject. As for measurement of the 

of estimating a number of dipoles in a short period of time. electromagnetic field distribution, an electroencephalograph 

Incidentally, the concrete content of method using the (i.e., EEG) is used with respect to measurement of electric 

pseudoinverse can be summarized as follows: 55 field distribution (i.e., brain waves) while a SQUID is used 

Singular value decomposition is performed on matrices with respect to measurement of magnetic field distribution, 

which are determined by position information of equivalent Normally, 20 to 150 measuring points are located on the 

current dipoles as well as arrangement of measuring points scalp of the subject. At the measurement, the system records 

of electromagnetic field distribution. Using information rep- coordinates of the measuring points as well as reference 

resenting a singular value, the method performs estimation 60 points of the head, for example, the nasion (i.e., the inter- 

of the equivalent current dipoles with respect their directions section of the internasal suture with the nasofrontal suture in 

and magnitude. the midsagittal plane) as well as the left and right ears, etc. 

There is a property that a square error at a time when the Incidentally, measurement of the shape of the head is 

aforementioned becomes minimal is equivalent to performed using the equipment of MRI (an abbreviation for 

magnitude of noise distribution. Using such a property, it is 65 "Magnetic Resonance Imaging") or X-ray CT scanner 

possible to improve noise resistance with regard to the (where "CT" is an abbreviation for "Computed 

dipole estimation method using the pseudoinverse. Con- Tomography"), for example. Instead of using an actual 
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shape of the head, it is possible to use an approximated 
model that approximates the head as a sphere, for example. 
Using such a sphere model may bring a small reduction in 
precision of estimation, however, there is an effect that load 
of calculations can be reduced. 5 

In step 2, the shape of the head of the subject is approxi- 
mated using polyhedrons (e.g., a set of triangles) on the basis 
of head shape data which are read into the system in step 1. 
Thus, the system creates grid points. FIG. 3 shows an 
example of a head model corresponding to an approximated 10 
shape of the head using a set of triangles. This head model 
is reconstructed from a set of MRI pictures of the head. As 
an active area model of the brain, a unit dipole having 
magnitude of 1.0 is located inside of each triangle. Herein, 
a number of model dipoles is represented by *m' (where 'm* 15 
is an integer arbitrarily selected). In this case, it is possible 
to perform calculations while treating a direction of the 
dipole as an unknown parameter. Instead, it is possible to 
calculate using dipoles with fixed directions perpendicular to 
the triangular surfaces. The above setting of the direction of 20 
the dipole is made based on the physiological knowledge 
that the axon of the neuron existing in the cortex is roughly 
perpendicular to the surface of the cortex. For this reason, 
the present embodiment is described under a precondition 
that the direction of the dipole is fixed perpendicular to the 25 
surface of the triangle. Under such a precondition, the 
present embodiment performs estimation to obtain the mag- 
nitude of the dipole only. Of course, it is possible to perform 
estimation in accordance with the same method of the 
present embodiment under a condition where the direction 30 
of the dipole is unknown. 

After completion of creation of the grid points in step 2, 
the system proceeds to step 3 so as to estimate dipole 
distribution by using the pseudoinverse. Incidentally, at a 
first circulation of the loop or in a state where the system 35 
does not meet the condition to get out from step 9 (which 
will be described later), it can be said that noise distribution 
has not been calculated yet. So, the aforementioned equation 
13 is used for estimation of the dipole distribution q. After 
completion of the estimation of the dipole distribution, the 40 
system calculates a square error Ek between measured 
electromagnetic field distribution of y and electromagnetic 
field distribution which is produced by the dipole in accor- 
dance with an equation as follows: 

45 

3=F*q [Equation 18] 

Thereafter, the system proceeds to step 4. In the sum of the 
equation 13 used for the estimation of the dipole distribution 
in the step 3, it is possible to estimate the dipole distribution 
in consideration of a singular value which is obtained by 50 
performing singular value decomposition of the matrix F. 
For example, the sum is taken with respect to the singular 
values X, such that Xi/X^xx, where Xq is a maximum singular 
value and a is a certain positive number. Thus, it is possible 
to avoid occurrence of divergence of the estimated values. If 55 
the estimation of the noise distribution is completed, the step 
3 uses the equation 17 to estimate the magnitude of the 
equivalent current dipole. In consideration of the estimated 
noise distribution, it is possible to reduce an effect of the 
noise by changing the equation 13 with the equation 17 with 60 
respect to the content of the sum. 

In step 4, a decision is made as to whether the noise 
distribution is estimated or not. In other words, the system 
performs checking as to whether the condition to get out 
from the loop is satisfied at step 9 or not. If the noise 65 
distribution is estimated, the system proceeds to step 12 to 
output the dipole distribution q as well as a number 'm' of 



the dipoles. Then, the system ends the proceeding. If the 
noise distribution is not estimated, the system proceeds to 
step 5. 

In step 5, the system performs calculations of the priority 
order to determine a place for each of elements of the dipole 
distribution q which is estimated by the step 3. So, the 
system sorts the elements based on their places of the 
priority order. Herein, the sorting is performed with respect 
to indexes of the elements of the dipole distribution q. For 
example, the elements are sorted as follows: 

where n is a permutation operator for the sorting of the 
elements. Accompanied with such an operation, ¥ £j - (i.e., an 
element ij of the matrix F) is permutated by F^. Thus, it 
is possible to execute the matrix calculations, required in the 
step 3, without changes at all. As the method of the sorting 
based on the priority order, there are provided two examples 
as follows: 

i) The absolute value of the magnitude of the dipole is 
used as the place of the priority order for the sorting. 
Thus, the elements are sorted into an order to arrange 
absolute values from the highest to the lowest. 

ii) A square error between measured electromagnetic field 
distribution and electromagnetic field distribution 
which is produced by a single dipole is used as the 
place of the priority order for the sorting. Thus, the 
elements are sorted into an order to arrange square 
errors from the lowest to the highest. 

After completion of the step 5, the system proceeds to step 
6. In the step 6, the system calculates an evaluation function 
based on the square error E k which is calculated in the 
step 4 as well as V which is a number of measuring points 
and 'm' which is a number of dipoles. As the evaluation 
function for example, it is possible to use the equations 
as follows: 

i) Use the structural risk defined by the aforementioned 
equation 16. 

ii) Use the description length L MDL defined by a first 
mathematical form of equation 20. 

iii) Use the Akaike information criterion defined by 
a second mathematical form of equation 20. 



nlogHT+- — logp 



[Equation 20] 



Laic = -IogE + m 



Incidentally, the above mathematical form defining the 
description length L MDL can be assumed as a simplified form 
of the equation 16. By using the description length h MI)L as 
the evaluation function it is possible to reduce an amount 
of calculations for the evaluation function. Incidentally, the 
theoretical basis for estimation of a number of parameters by 
using the description length is provided by "paper 5" cor- 
responding to Rissanen, "Modeling by Shortest Data 
Description", Automatica, Vol. 14, pp. 465^71, 1978, for 
example. 

In addition, the theoretical basis for estimation of a 
number of parameters by using the Akaike information 
criterion shown by the equation 20 is provided by 
"paper 6" corresponding to a book written by Mr. Ishiguro 
and his members which is entitled "Information Theoretical 
Statistics" published by Kyoritsu Publication Co. of Japan 
on 1983, for example. 

After completion of the step 6, the system proceeds to step 
7 wherein a decision is made as to whether the loop counter 
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variable k is zero or not. If the variable k is zero, the system 
proceeds to step 8. If not, the system proceeds to step 9. In 
step 8, '1* is added to the loop counter variable k. Then, the 
system proceeds back to the aforementioned step 3. In step 
9, the system compares the "present" evaluation function L k 
with a previous evaluation function L^, wherein is 
produced at present execution of the loop while L^, is 
produced at previous execution of the loop. So, a decision is 
made as to whether an inequality of L Jt . 1 >L Jt is established or 



dipoles is three. Of course, this invention is capable of 
estimating a number of dipoles even if a number of test 
dipoles is one, two or four. 

We have made another study in availability of this inven- 
tion with respect to the case of the distributed active area of 
the brain while the noise is contained in measurement data. 
Herein, estimation is performed using the pseudoinverse. 
Simulation data subjected to the estimation are produced 
using a "distributed" dipole model of FIG. 4. In addition, the 



not. There is a tendency that a value of the evaluation 10 random noise generated on the basis of the Normal distri 



function becomes small as the estimation is appropriate. So, 
the inequality of L Jt . 1 >L A indicates that calculations of the 
estimation come to convergence smoothly. For this reason, 
under the inequality of L^. 1 >L Jfc) the system proceeds to step 
10 so as to perform estimation of the dipole distribution 
while changing a number of dipoles. On the other hand, an 
inequality of L*.., ^L* indicates that a result of the estimation 
is obtained for the time being. So, the system proceeds to 
step 11 to obtain a final result of the estimation in consid- 
eration of the effect of the noise. 

In step 10, a number of dipoles is changed on the basis of 
the priority order which is determined by the aforemen- 
tioned step 6. As the method of changing the number of the 
dipoles, there are provided two methods as follows: 



15 



20 



bution is added to the simulation data, and the magnitude of 
the noise is set to 20% of the simulation data. Results of the 
estimation are shown in FIG. 5A and FIG. 5B. Specifically, 
FIG. 5A shows results of the estimation based on the method 
of this invention while FIG. SB shows results of the esti- 
mation truly based on the conventional method using the 
pseudoinverse. In FIG. 4 and FIG. 5A, intensity of the 
equivalent current dipole (see vertical axis) is represented in 
an order of xlO -1 Angstrom-meter [A°m]. On the other 
hand, the intensity of the current dipole of FIG. 5B is 
represented in an A°m unit. So, the intensity of the equiva- 
lent current dipole of FIG. 5A differs from that of FIG. SB 
by 10 10 times or so. This invention is capable of detecting a 
center location of the original dipole distribution. In 



of q^ having a smallest value is deleted, so that *r 
which is a number of dipoles is decreased by '1*. 



i) Among elements of the dipole distribution q, an element 25 contrast, the conventional methods cannot detect a center 

*m' location of the dipole distribution because the dipole distri- 
bution is affected by the noise and its estimation is subjected 

ii) Among elements of the dipole distribution q, the to divergence. 

system searches an element having a highest place of As this invention may be embodied in several forms 

the priority order within a group of dipoles which are 30 without departing from the spirit of essential characteristics 

not used until now. This element is added to the dipole thereof, the present embodiment is therefore illustrative and 

distribution, so that 'm* which is a number of dipoles is not restrictive, since the scope of the invention is defined by 

increased by ' 1 \ the appended claims rather than by the description preceding 

After changing the number of the dipoles, the system them, and all changes that fall within metes and bounds of 

proceeds to step 8 wherein l V is added to the loop counter 35 the claims, or equivalence of such metes and bounds are 

variable k. Thereafter, the system proceeds back to the therefore intended to be embraced by the claims, 

aforementioned step 3. What is claimed is: 

In step 11, when estimation is performed using a certain 1. An electrophysiological activity estimation method 

number of dipoles which brings L k to coincide with a which inputs data of measured electromagnetic field distri- 

minimum value, the square error is used as an estimated 40 bution measured on a surface of a living body and data of a 

value of the magnitude of the noise distribution of ||z||. After shape of the living body so as to assume at least one 

completion of the step 11, the system proceeds back to the equivalent current dipole in the living body as a source of the 

step 3. measured electromagnetic field distribution, said electro- 

By executing the aforementioned steps, the system pro- physiological activity estimation method comprising the 

ceeds to step 12 at last, wherein the system outputs results 45 steps of: 



of the estimation, Le., the dipole distribution (i.e., directions 
and magnitude of dipoles) and a number of dipoles. Then, 
the system ends processing of the flowchart of FIG. 1. 

Next, a description will be given with respect to results of 
evaluation for availability of the present invention used to 50 
perform estimation of a number of dipoles. Herein, simula- 
tion data used for the above evaluation is produced using 
three test dipoles. So, the method of this invention is used to 
perform the estimation while changing a number of model 
dipole(s) from one to eleven. Results of the evaluation are 55 
given by variations of an evaluation function L, which are 
shown in a graph of FIG. 3 wherein a vertical axis represents 
a value of the evaluation function L while a horizontal axis 
represents a number of model dipoles. Incidentally, the 
aforementioned structural risk is used as the evaluation 60 
function L. FIG. 3 shows that the value of the evaluation 
function L becomes minimal when a number of the test 
dipoles coincides with a number of the model dipoles. 
Namely, it is understood from the content of FIG. 3 that the 
new evaluation criterion (e.g., structural risk) employed by 65 
this invention is effective in estimation of a number of 
dipoles. FIG. 3 shows the case where a number of test 



performing first estimation with respect to a direction and 
magnitude of the equivalent current dipole such that an 
error between the measured electromagnetic field dis- 
tribution and a theoretical value of electromagnetic 
field distribution caused by the equivalent current 
dipole becomes minimal; 

detennining a priority order for the estimated equivalent 
current dipole; 

calculating an evaluation function on the basis of a 
number of measuring points used for measurement of 
the measured electromagnetic field distribution and a 
number of equivalent current dipoles; 

changing the number of equivalent current dipoles based 
on the evaluation function; 

performing second estimation with respect to noise dis- 
tribution on the basis of the evaluation function; and 

outputting the direction and magnitude of the equivalent 
current dipole as well as the number of equivalent 
current dipoles. 

2. An electrophysiological activity estimation method as 
defined in claim 1 wherein a structural risk calculated by the 



04/22/2004, EAST Version: 1.4.1 



6,073,1 

13 

error, the number of measuring points and the number of 
equivalent current dipoles is used as the evaluation function. 

3. An electrophysiological activity estimation method as 
defined in claim 1 wherein a description length calculated by 
the error, the number of measuring points and the number of 
equivalent current dipoles is used as the evaluation function. 

4. An electrophysiological activity estimation method as 
defined in claim 1 wherein Akaike information criterion 
calculated by the error, the number of measuring points and 
the number of equivalent current dipoles is used as the 
evaluation function. 

5. An electrophysiological activity estimation method as 
defined in claim 1 further comprising the step of: 

performing deletion of equivalent current dipoles from a 
lowest place of the priority order. 

6. An electrophysiological activity estimation method as 
defined in claim 1 further comprising th e step of: 

performing addition of equivalent current dipoles from a 
highest place of the priority order. 

7. An electrophysiological activity estimation method as 
defined in claim 1 wherein the priority order is determined 
using an absolute value of the equivalent current dipole. 

8. An electrophysiological activity estimation method as 
defined in claim 1 wherein the priority order is determined 
using the error between the measured electromagnetic field 
distribution and the theoretical value of the electromagnetic 
field distribution caused by the equivalent current dipole. 

9. An electrophysiological activity estimation method as 
defined in claim 1 wherein after the second estimation of the 
noise distribution based on the evaluation function, the first 
estimation is performed again in such a way that a matrix 
determined by position information of the equivalent current 
dipole and arrangement of the measuring points of the 
electromagnetic field distribution is subjected to singular 
value decomposition, so that the direction and the magnitude 
of the equivalent current dipole are estimated using the noise 
distribution as well as information of the error. 

10. An electrophysiological activity estimation method as 
defined in claim 1 wherein the first estimation is performed 
in such a way that a matrix determined by position infor- 
mation of the equivalent current dipole and arrangement of 
the measuring points of the electromagnetic field distribu- 
tion is subjected to singular value decomposition, so that the 
direction and the magnitude of the equivalent current dipole 
is estimated using information of a singular value. 

11. An electrophysiological activity estimation method 
comprising the steps of: 

inputting data representing measured electromagnetic 
field distribution which is measured at measuring 
points on a surface of a selected part of a living body 
as well as data representing a shape of the selected part 
of the living body; 

creating grid points with respect to the surface of the 
selected part of the living body; 

performing estimation of dipole distribution for equiva- 
lent current dipoles each of which is assumed as a 
source of the electromagnetic field distribution within 
the selected part of the living body; 

determining places of a priority order with respect to the 
equivalent current dipoles respectively; 

sorting the equivalent current dipoles in accordance with 
the priority order; 



14 

performing calculation to produce an evaluation function 
in accordance with a square error between the mea- 
sured electromagnetic field distribution and a theoreti- 
cal value of electromagnetic field distribution which is 
caused to occur due to the equivalent current dipole as 
well as a number of the measuring points and a number 
of the equivalent current dipoles; and 
changing the number of the equivalent current dipoles in 
accordance with the priority order under a condition 
where a result of the calculation comes to convergence, 
so that the aforementioned steps are repeated while 
changing the number of the equivalent current dipoles. 

12. An electrophysiological activity estimation method as 
15 defined in claim 11 wherein the estimation is performed 

using a pseudoinverse while a structural risk is used as the 
evaluation function. 

13. An electrophysiological activity estimation method as 
defined in claim U further comprising the step of: 

20 performing estimation of noise distribution if the result of 
the calculation departs from the convergence, so that 
the dipole distribution is estimated using the noise 
distribution. 

25 14, An electrophysiological activity estimation method as 
defined in claim 11 further comprising the step of: 

performing estimation of noise distribution corresponding 
to a square error calculated with respect to a number of 
equivalent current dipoles which gives a m inimum 
30 value to the evaluation function, so that the dipole 
distribution is estimated using the noise distribution. 

15. An electrophysiological activity estimation method as 
defined in claim 11 wherein the dipole distribution is rep- 
resented by a direction and magnitude of the equivalent 

35 current dipole as well as the number of the equivalent 
current dipoles which are estimated in such a way that the 
square error becomes niinimal. 

16. An electrophysiological activity estimation method as 
defined in claim U wherein the dipole distribution is rep- 

40 resented by a direction and magnitude of the equivalent 
current dipole which are estimated using a singular value 
which are obtained by subjecting a matrix determined by a 
position of the equivalent current dipole and arrangement of 
the measuring points to singular value decomposition. 
45 17. An electrophysiological activity estimation method as 
defined in claim U wherein the number of the equivalent 
current dipoles is changed in such a way that an equivalent 
current dipole having a smallest value within the equivalent 
current dipoles of the dipole distribution is deleted, so the 
number of the equivalent current dipoles is decreased by 
one. 

18. An electrophysiological activity estimation method as 
defined in claim U wherein the number of the equivalent 
current dipoles is changed in such a way that an equivalent 

55 current dipole having a highest place of the priority order 
within equivalent current dipoles which are not used for the 
estimation is added to the dipole distribution, so the number 
of the equivalent current dipoles is increased by one. 

19. An electrophysiological activity estimation method as 
60 defined in claim 11 wherein the selected part of the living 

body corresponds to a head of a human. 

* * * * * 
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